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Abstract 

An optimal estimation inverse method is presented which can be used to retrieve simultaneously 
vertical profiles of temperature and specific humidity, in addition to surface pressure, from 
satellite-to-satellite radio occultation observations of the Earth's atmosphere. The method is a 
non-linear, maximum a posteriori technique which can accommodate most aspects of the real 
radio occultation problem and is found to be stable and to converge rapidly in most cases. The 
optimal estimation inverse method has two distinct advantages over the analytic inverse method 
in that it accounts for some of the effects of horizontal gradients and is able to retrieve optimally 
temperature and humidity simultaneously from the observations. It is also able to account for 
observation noise and other sources of error. Combined, these advantages ensure a realistic 
retrieval of atmospheric quantities. 

A complete error analysis emerges naturally from the optimal estimation theory, allowing a full 
characterisation of the solution. Using this analysis a quality control scheme is implemented 
which allows anomalous retrieval conditions to be recognised and removed, thus preventing gross 
retrieval errors. 

The inverse method presented in this paper has been implemented for bending angle 
measurements derived from GPS/MET radio occultation observations of the Earth. Preliminary 
results from simulated data suggest that these observations have the potential to improve NWP 
model analyses significantly throughout their vertical range. 
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1. Introduction 

Radio occultation (RO) experiments have played 
a prominent role in the NASA programme for so- 
lar system exploration for more than two decades 
and have contributed to studies of the atmosphere 
of Mars [Fjeldbo and Eshleman, 1968], Venus [Fjeldbo 
and Kliore, 1971], Jupiter [Kliore et al, 1975], Saturn 
[Lindal et al, 1985], Uranus [Lindal et al, 1987], and 
Neptune [Lindal 1992]. This method of radio occulta- 
tion uses a receiver on Earth and a satellite occulted 
by a planetary atmosphere (which may occur from 
either a fly- by or by a satellite orbit of the planet). 
Suitably accurate atmospheric RO measurements of 
the Earth's atmosphere became possible with the ad- 
vent of the Global Positioning System (GPS), but it 
was not until the late 1980s— early 1990s that the po- 
tential of RO using the GPS was widely appreciated 
(e.g. Gurvich and Krasil'nikova, [1990]). 

The radio occultation method used to sound the 
Earth's atmosphere is different from that used by the 
planetary experiments, in that both the receiver and 
the transmitters are orbiting the planet. 

Data from the prototype GPS space-borne re- 
ceiver, GPS/MET, launched in April 1995, confirmed 
the potential of obtaining accurate, global observa- 
tions of the Earth's atmosphere from the radio occul- 
tation technique. Temperature comparisons between 
early results from the GPS/MET receiver and collo- 
cated radiosondes and numerical weather prediction 
(NWP) model analyses showed good agreement (e.g. 
Kursinski et al, [1996]; Ware et al, [1996]). 

The analytic method of inverting radio occulta- 
tion measurements to obtain meteorological param- 
eters (i.e. the method used to sound other planetary 
atmospheres) involves the use of an integral trans- 
form, using the assumption of a horizontally homo- 
geneous atmosphere, to obtain a profile of refractiv- 
ity (as a function of geometric height) [Fjeldbo and 
Eshleman, 1968]. The hydrostatic relation is used 
to obtain pressure and temperature from refractiv- 
ity via density. For the Earth's atmosphere, where 
a reasonable prior knowledge of horizontal gradients 
is available, the analytic inversion does not represent 
the most suitable method since inadequate modelling 
of such gradients can cause large retrieval errors (e.g. 
Ahmad and Tyler, [1998]). 

Eyre, [1994] addresses this issue and suggests a sta- 
tistically optimal retrieval approach, using variational 
methods, to enable the direct assimilation of bending 
angle or refractivity (Healy and Eyre, [1999] investi- 



gate the latter quantity). Zou et al, [1995] also looked 
at the impact of atmospheric radio refractivity mea- 
surements using a 4-D variational data assimilation 
approach. Their results showed that the measure- 
ments were effective in recovering the vertical pro- 
files of water vapour, and found that the accuracy of 
the derived water vapour field was significantly bet- 
ter than that obtained through the analytic retrieval 
technique. The assimilation of these measurements 
were also shown to provide useful temperature infor- 
mation. There have also been several numerical ex- 
periments which have assessed simulated GPS/MET 
refractivity measurements to predict cyclonic distur- 
bances (e.g. Kuo et al, [1998]), and have concluded 
that these measurements are likely to have a signif- 
icant impact on short-range operational NWP, with 
the caveat that the number of GPS receivers will have 
to be increased before the full potential impact of this 
measurement could be realised. 

In this paper we utilise a non-linear optimal esti- 
mation technique which is implemented and validated 
using an ensemble of simulated retrieval scenarios, us- 
ing the bending angle quantity as the 'observation'. 

Section || outlines the details of the RO method- 
ology for the Earth's atmosphere necessary to derive 
the bending angle quantity, recalls the analytic in- 
verse method and discusses the impetus for pursuing 
an alternative inverse method. In section |^ we out- 
line the theory for the non-linear optimal estimation 
inverse method and give details of its implementation 
for GPS RO observations. Section ^ is devoted to 
details of the validation of the optimal estimation in- 
verse method with reference to GPS RO observations. 
The sensitivity of the inverse model assumptions is 
also investigated, and we conclude the paper with a 
discussion of the results obtained. 

2. Radio occultation measurements of 
Earth's atmosphere 

Kursinski et al, [1997] give a detailed description 
of the method used to measure the RO atmospheric 
obscrvables; the following section gives a summary of 
the theory, assuming no external encryption of the 
signals. 

The GPS satellites transmit on two L-band radio 
frequencies^. Assuming a continuous link between the 
receiver and transmitter, when the receiver passes be- 
hind the atmosphere with respect to a GPS transmit- 

iNamely, LI: 1.57542 GHz and L2: 1.2272 GHz. 
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ter the signal travels through the atmosphere and is 
refracted in response to variations of refractive index 
along its path. This refraction causes the ray to travel 
over a longer path than it would in the absence of the 
atmosphere, in accordance with Fermat's principle of 
least time, which subsequently causes an atmospheric 
time delay in the received signal. 

The Doppler shift of the signal is calculated from 
the additional atmospheric delay (the derivative of 
the phase delay). Using the geometry and notation of 
figure |l|, the Doppler shift /d of the carrier frequency 
fa measured by the receiver is given by: 



fd — fo 



(VT ■ riT + Vr • Hr) 



(1) 



where vt,r are transmitter and receiver velocity vec- 
tors, riTjR are path direction vectors of the transmit- 
ter and receiver, and c is the velocity of light. 

There are also relativistic terms which need to 
be considered in equation |l| (due to different grav- 
itational potentials and higher order corrections for 
spacecraft velocity) but these can be eliminated based 
on knowledge of orbital geometry and the Earth's 
gravity field [Kursinski et al, 1997]. Note that the rel- 
ative positions and velocities of the two satellites can 
be calculated very accurately using available tracking 
data, which is independent of radio occultation data. 

By specifying radial and tangential components of 
the velocity of satellite i in the plane coinciding with 
the ray trajectory as and v*, and taking into ac- 
count Snell's law, the angles 0r and 0t can be calcu- 
lated from the following relations: 



/d 



foe ^(VrCOS( 



■ Vqn COS 0T 



sin</)R - sin(/)T) 



(2) 



The cumulative effect of the atmosphere on the ray 
path can be expressed in terms of the total refractive 
bending angle, e, as a function of the impact param- 
eter, a. The impact parameter may be defined as 
the perpendicular distance between the local curva- 
ture of the Earth at the tangent point of the ray and 
the asymptotic straight line followed by the ray as it 
approaches the atmosphere. 

From Bouguer's rule [Born and Wolf, 1993], and 
the geometry defined by figure ^ e(a) can be calcu- 
lated thus 



tr sm 0R = tt sm (p^ = a 

e — <j)R + 4>T + 6- TT 



(3) 
(4) 



It is this rule that introduces the assumption of spher- 
ical symmetry (nr sin (f) = a, where a is a constant 
along the ray path), i.e. a horizontally homogeneous 
atmosphere. Departures from this assumption can 
introduce significant errors if not properly accounted 
for. These errors have been studied by Ahmad and 
Tyler, [1999] and Healy, [1999], but are not addressed 
in the study presented here. 

Measurements of the time delay become possible 
for neutral atmospheric sounding when the GPS sig- 
nal begins to transect the mesosphere at an altitude of 
about 85 km; at this altitude the atmospheric phase 
delay is about 1 mm (3 x 10^^^ s) which can be ob- 
served by the LEO GPS receiver [Ware et al, 1996]. 
Further information about the measurement charac- 
teristics may be obtained from Kursinski et al, [1997]. 

2.1. The analytic inverse method 

Using an Abel integral transform (equation 
or making use of a similar integral transform when 
applying Fresnel diffraction theory [Mortensen and 
H0eg, 1998], these bending angle measurements can 
be inverted to obtain a profile of refractivity. For com- 
pleteness sake, the 'forward' Abel integral transform 
(equation ||) is presented alongside the 'inverse' Abel 
integral transform: 



Inn(x) = 



1 



+ 00 



e{a){a^ 



'r'^^da (5) 



e(a) = -2a^ °°(^^)(.^-ar^/^dx, (6) 

where n is the refractive index and x is the refractive 
radius (i.e. x = rn). 

Geometric height levels, z, can be obtained from 
the refractive index profile, as a function of impact 
parameter, and the local radius of curvature, Rc, thus 



X 

Z = i?c 

n 



(7) 



Because refractive index is near unity, refractivity 
N is used to describe the refractive medium which is 
given by iV = (n — 1) X 10^. 

Refractivity is affected primarily by air density 
(dependent on pressure and temperature) and water 
vapour density, thus the measurement contains infor- 
mation about both. Equation ^ describes this rela- 
tionship. 
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N = 



ClPa . C2Pw 



(8) 



where Pa is the total atmospheric pressure (hPa), P^f 
is the partial pressure of water vapour, T is the tem- 
perature (K), and ci and C2 represent constants of 
proportionality, whose values are 77.6 (KhPa^^) and 
3.73 X lO'^ (K^hPa-i), respectivelv. The form of the 
dry and moist terms in equation ^ is from Smith and 
Weintraub, [1953]. 

Refractivity is also affected by charged particles in 
the ionosphere and the scattering by water droplets 
suspended in the atmosphere. The first-order iono- 
spheric contribution to refractivity can be removed 
by combining the two GPS signals (described by 
Vorob'ev and Krasil'nikova, [1997]), leaving higher- 
order terms, and the scattering contribution is found 
to be negligible compared to the contribution due to 
air and water vapour density [Kursinski et al, 1997]. 

There is no measurement information to allow the 
separation of the effects of temperature and water 
vapour, and therefore these quantities can be re- 
trieved only using prior information. If the abso- 
lute humidity is judged to be small (e.g. in the cold- 
est regions of the troposphere and stratosphere, with 
temperatures less than 250 K), it may be neglected, 
and density calculated from refractivity. The hydro- 
static relation can be used to calculate values of pres- 
sure, and hence temperature. However, if humidity 
is judged to be significant, then an iterative process 
may be used to calculate temperature/humidity if an 
a priori profile of humidity/temperature is used. This 
prior information can be taken from various sources 
such as collocated NWP model output. 

The inability of this inverse method to account for 
horizontal refractivity inhomogeneities, and the sub- 
optimal way this method retrieves temperature and 
humidity with prior values, represents two key disad- 
vantages of this method, and form part of the impetus 
to develop a new inverse method. 

The hydrostatic relation, used to compute values 
of pressure on the retrieved height levels, requires 
an assumed pressure value at a particular geometric 
height level. A variety of methods have been imple- 
mented to tackle this initial value problem. Kursin- 
ski et al, [1996] assumed a temperature of 260 K at 
50 km. This method has the problem that if the as- 
sumed temperature is inconsistent with the measure- 
ments an error is introduced, which decreases expo- 
nentially with depth. More elaborate methods (e.g. 
Rocken et al, [1997] and Steiner et al, [1999]) initialise 



the GPS/MET retrieval at some high altitude (e.g. 
100 km) using climate model data, and combine the 
measurements and model data to minimise downward 
propagation of errors. In principle, the method pre- 
sented in this paper also combines the model data 
and the observations but achieves this is an optimal 
way. The method presented also has the advantage 
of a straightforward error characterisation. 

Since GPS RO observations sometimes reach near- 
surface altitudes (i.e. less than 1 km from the sur- 
face), surface pressure is also retrieved using the op- 
timal estimation inverse method. 

3. The optimal estimation inverse 
method 

The method outlined here is also known as one- 
dimensional variational data analysis. 

The main advantage of this method is that it pro- 
vides simultaneous estimates of temperature and hu- 
midity profiles that are statistically optimal, given 
prior estimates from an NWP model (together with 
their error covariances). It also provides a framework 
for assessing the error characteristics of the estimates. 

In this study only the 1-d problem has been stud- 
ied. However, the errors introduced by the neglect of 
the horizontal gradients have been estimated and al- 
lowed for as part of the error budget (see section 3.2 
"Forward modelling errors" ) . 

3.1. Theory 

A brief description of the theory used in optimal 
estimation in presented here; a more detailed descrip- 
tion may be found in Rodgers, [1976] and Rodgers, 
[1990]. 

For brevity, the observation noise, the error associ- 
ated with any forward modelling parameters and the 
forward model error (which includes the representa- 
tiveness error [Lorenc, 1986] ) will be accounted for in 
one vector which will be denoted by e and its ensem- 
ble characteristics described by the covariance matrix 
E. 

The rationale behind optimal estimation is to min- 
imise a cost functional J(x) (or to solve VxJ(x)=0), 
which measures the degree of fit of estimates of the 
atmospheric state to the measurements and to some 
prior information, and possibly to some other physical 
or dynamical constraints. In this case J(x) is given 

by 
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J(x) = (y°-y(x)fE-i(y°-y(x)) + 

(x-x'^)Tc-i(x-x^^) (9) 

where x*^ and x represent the background and up- 
dated state vectors, respectively; y° and y(x) repre- 
sent the observation vector and the estimated obser- 
vation vector calculated from the state vector {Eyre, 
[1994]), respectively; and C represents the back- 
ground error covariance matrix. 

There are a number of methods available to min- 
imise J(x): the scheme described here uses the Leven- 
berg-Marquardt iterative method [e.g. Press et al, 
1992]: 

Xi+i = x'^-t-((l+7)C-i+KTE-iK)-i 
[(KTE-i(y°-y(xO)) + 
(7C-i + KTE-iK)(xi-x^)] (10) 

where K is Vxiy(xi), 7 is a non-dimensional weight- 
ing factor 0, and all other variables are as before. 

Using the optimal estimation theory it is possible 
to obtain an error covariance for the retrieved prod- 
ucts. Indeed, it can be argued that the retrieved prod- 
ucts are of limited value without an estimate of their 
uncertainty. The solution error covariance S is given 
approximately (i.e. at the linear limit) by 

S = (C-i +K'rE-i K)-^ (11) 

The solution error covariance can then be com- 
pared with the prior error covariance to ascertain how 
the retrieval has improved upon the prior knowledge 
of the atmospheric state. 

3.2. Implementation of the optimal 
estimation inverse model 

This subsection describes in detail the components 
of equation ^. 

The background state vector and its uncer- 
tainty covariance matrix In this case the back- 
ground knowledge of the atmosphere state x*^ was 
obtained from short-range forecasts provided by the 
UKMO unified model [Cullen, 1993]. The model 
from which the data are derived had 19 levels, which 
were expressed on hybrid-sigma pressure coordinates 

■^for increasing values of 7 this minimisation method degen- 
erates into the method of steepest descent. 



(surface— 10 hPa). The global model had a resolution 
of 0.833° (I8O7217) latitude and 1.25° (360°/288) 
longitude. 

The data used are 6-hour forecasts which have been 
interpolated to occultation event positions, using the 
mean latitude and longitude of each occultation^. 

These 19 levels are linearly interpolated (in In- 
pressure) onto the state vector levels used for TOVS 
retrievals [Eyre, 1990]. CIRA chmatology [CIRA, 
1986] is assumed above the UKMO model account- 
ing for the latitudinal and seasonal variation of the 
profile. This climatology provides a reasonable prior 
and first guess information in the upper stratosphere. 

For the forecast error covariance matrix C (de- 
scribed in Eyre, [1989]), lower atmospheric values 
(surface— 50 hPa) were generated from radiosonde— forecast 
difference statistics and upper stratosphere values 
were found by regression from the levels provided 
[Eyre, 1989]. 

The radio occultation retrieval uses 40 tempera- 
ture elements, 15 ln(specific humidity) elements^ and 
a surface pressure element from this forecast error co- 
variance matrix. The temperature and ln(specific hu- 
midity) inter-quantity covariance values have been set 
to zero. These inter-quantity covariances are not well 
known and assuming zero covariance between them is 
more conservative than an erroneous covariance. The 
surface pressure element is uncorrelated with both 
temperature and ln(specific humidity). 

Since CIRA climatology is used to form the a pri- 
ori (and the first-guess) it is necessary to consider 
the errors that may be attached to such information. 
In general, if the standard deviation values from the 
diagonal elements of the forecast error covariance ma- 
trix are smaller than the uncertainties assumed for the 
climatology, then the climatological errors are used at 
the levels in the upper atmosphere described by the 
CIRA climatology (off-diagonal elements remain the 
same): at latitude 6, for \9\ > 45° a— 15 K (winter) 
and (7=5 K (summer); and for \9\ < 45° (7=5 K. 

The values for the diagonal elements of the UKMO 
forecast error covariance matrix are shown by figure 

I 

As expected for temperature, the lower atmosphere 

^The mean latitude and longitude of an occultation corre- 
sponds typically to altitudes in the lower stratosphere/upper 
troposphere. 

^Specific humidity is expressed as the natural logarithm of 
specific humidity since forecast errors in this quantity are more 
constant than in specific humidity. 
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forecast errors are reasonably small (of the order of 
1.5 K) and increase as a function of altitude. The 
actual values have been developed over recent years 
at the UKMO and reflect the average error in six-hour 
forecasts. 

The observation vector and its error covari- 
ance matrix The observation vector y° contains 
bending angle measurements as a function of impact 
parameter (section ||). 

In practice, atmospheric phase delay measurements 
from the GPS/MET receiver are low-pass filtered to 
reduce noise. The cut-off frequency of the filter is 
tuned to pass phase variations corresponding to verti- 
cal scales of 2 to 3 km in the stratosphere and approx- 
imately 200 m in the lower troposphere [Rocken et al, 
1997]. From the phase observable, Doppler shifts (and 
subsequently bending angle profiles) for the two GPS 
signals are computed. First-order ionospheric effects 
are removed from the data by combining the two sig- 
nals to form a single corrected profile [Vorob'ev and 
KrasiVnikova, 1994]. Typically, after filtering, there 
are 100—200 neutral atmosphere bending angle mea- 
surements, which span a vertical range of approxi- 
mately 0.5—60 km. 

In addition to the observation error covariance ma- 
trix consisting of observation noise estimates, errors 
from the forward modelling and forward model pa- 
rameters are considered. 

Observation noise Observation errors are cre- 
ated by the hardware of the measurement system and 
by the pre-processing of the observations. The ob- 
servation noise estimates are taken from the work 
described in Luntama, [1997]. They include ther- 
mal noise, residual errors from the ionospheric cor- 
rection, local multipath (distortions when the trans- 
mitted signal is refiected from a surface near the sig- 
nal propagation path), orbit determination accuracy, 
and clock instabilities of low-Earth-orbit receiver and 
GPS satellites and ground stations. 

These observation error estimates are based on 
phase noise levels during the measurement or esti- 
mated from other noise sources during a radio occul- 
tation event. The effect from satellite clock errors 
and from the selective availability military encryp- 
tion process were found to be negligible by assum- 
ing a differencing decryption technique (involving the 
differencing of the several sets of signals, using the 
satellites and the ground stations) [Kursinski et al, 
1997]. The only exception is the residual error from 
the ionospheric correction which was obtained from 
refractivity error estimates published in Kursinski et 



al, [1997] and mapped to bending angle space using a 
forward model [Luntama, 1997]. 

The observation error estimate^ used in the op- 
timal estimation technique are shown by figure |[ 
These estimates represent normal atmospheric condi- 
tions with a relatively small multipath error (3 mm) 
and normal ionospheric conditions. Panel (a) shows 
that there are a number of noise contributions of com- 
parable size in the lower atmosphere. At altitudes 
above 30 km the residual ionospheric correction er- 
ror begins to dominate the total bending angle error 
curve (panel (b)). 

These noise estimates are assumed to be fully in- 
dependent, i.e. their inter-level (and inter-quantity) 
covariance is zero. However, the bending angle mea- 
surements do contain a small, local correlation be- 
tween successive levels due to filtering of the phase 
measurements. Because this correlation is small, the 
diagonal form is a good approximation to the full ma- 
trix [J. P. Luntama: personal communication]. 

Real observations from the GPS/MET data used 
have been found to be noisier than theoretical esti- 
mates [Luntama, 1997]. Bending angle fluctuations 
in the upper stratosphere (of the order of 10^^ radi- 
ans) are present and are thought to be due to residual 
errors from the LEO satellite clock calibration in the 
differencing decryption technique (see "Observation 
noise") [Syndergaard, 1999]. As such, a suitable error 
is attached to reflect the upper atmosphere measure- 
ments. 

Forward modelling errors In this work the 
forward model used to map from state space to obser- 
vation space is described in Eyre, [1994] but applied 
to an atmosphere approximated as spherical symmet- 
ric about the given profile at the tangent point. Es- 
sentially the geophysical parameters are converted to 
refractivity as a function of height, and subsequently 
impact parameter using the local radius of curvature. 
The resulting profile is mapped into observation space 
using the 'forward' Abel integral transform (equation 

D- 

The two main forward modelling errors are due to 
the assumption of a horizontally homogeneous atmo- 
sphere and a representativeness error. 

Estimates for the first of these errors are obtained 
using a version of the forward model which can ac- 

^ These estimates have been computed by a nominal bend- 
ing angle profile defined using an exponential curve with a 
scale height of approximately 7 km [J. P. Luntama: personal 
communication] . 
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count for horizontal inhomogeneities in the plane of 
the ray path (described in Eyre, [1994]) and com- 
paring observation vectors with the version of the 
forward model which assumes a horizontal homoge- 
neous atmosphere [Palmer, [1998]]. Mid-latitude two- 
dimensional NWP fields (0—360°) were used to sim- 
ulate typical horizontal gradients in temperature and 
humidity. By considering small sections of the field 
at a time (typical of the horizontal resolution of ra- 
dio occultation measurements which is of the order of 
300 km), the complete field was traversed. Comput- 
ing the ensemble mean from the difference between 
the two versions of the forward models allowed a rea- 
sonable estimate of the forward modelling error in- 
curred by the spherical symmetry assumption to be 
computed. It is noted that this error estimate does 
not represent the true error in observation space since 
the forward model does not simulate the full error 
characterisation. Both Ahmad and Tyler, [1999] and 
Healy, [1999] consider bending angle errors from hori- 
zontal gradients for specific cases. However there is no 
published material that quantifies this error statisti- 
cally. Simulation with a full 3-d raytracer through the 
UKMO mesoscale model fields suggest that the errors 
are approximately 3% for ray paths near the surface, 
which is consistent with the value used in this work. 

An error arises from representing an intrinsically 
high resolution problem with a crude resolution. This 
type of error is often called a representativeness error 
and describes the error from the inability of NWP 
model vertical grids to represent small-scale atmo- 
spheric structure, which are evident in GPS/MET RO 
measurements [Kursinski et al, [1997]]. The method 
used to estimate this quantity is described by Healy, 
[1998], and is found to be of the order of 2% of the 
bending angle measurement in the troposphere and 
upper stratosphere, decreasing slightly in the middle 
stratosphere. This variation in the error is associated 
with the temperature variations in these region. 

The forward model used is based on geometric 
optics therefore does not account for atmospheric 
diffraction. However, Kursinski et al, [1997] have 
shown that the geometric optics assumption is suc- 
cessful in describing propagation characteristics above 
a certain diffraction limit, and diffraction is therefore 
not considered here. 

Forward model parameter errors Uncertain- 
ties associated with the physical constants used to 
model the physical system also cause modelling er- 
rors. The major physical constants used in the for- 
ward model are the refractivity coefficients (ci and C2 



in equation |^) and the local radius of curvature. 

The uncertainty of the refractivity coefficients do 
not represent the error associated with their values 
but the uncertainty of the measured quantity; for this 
reason the information will not be included in the to- 
tal observation error budget since it will result in a 
bias in the retrieval The local radius of curvature 
assumed for zero altitude is a parameter that is used 
to compute the bending angle observations from at- 
mospheric phase delay, therefore any error associated 
with this parameter will be present in the bending an- 
gle observations. This parameter is also used to com- 
pute geometric height levels from impact parameter 
levels. The uncertainty of this value is estimated to 
be approximately 100 metres [Kursinski et al, [1997]]. 

Total observation error The total observation 
error covariance E is constructed by adding the co- 
variance matrices from observation noise, forward 
modelling and forward model parameters. 

The diagonal terms have also been constrained not 
to fall below a minimum value to account for the noisy 
upper stratosphere measurements. 

Figure ^ shows how the standard deviation values 
of the principal diagonal from each error contribute 
to the observation error covariance matrix. 

The dominant source of error for the majority of 
the vertical range considered is the forward modelling 
error, i.e. representativeness error and horizontal in- 
homogeneity error, with the upper stratospheric noise 
limit providing the second largest contribution to the 
total error. In the upper stratosphere the total bend- 
ing angle reverts to the upper level minimum noise 
used. At near-surface altitudes, the forward model 
parameter error, i.e. local radius of curvature, is sig- 
nificant. The overall effect from the local radius of 
curvature decreases exponentially due to the hydro- 
static relation. 

3.3. Convergence and quality control 

The method used to judge convergence relies upon 
values of the cost function, i.e. if the relative change is 
smaller than a specified value (0.5%) then the solution 
is determined to have converged. This method alone 
is found to be a good indicator of convergence in this 
case. 

For the work presented in this paper, the maxi- 
mum number of iterations considered is 10; if the so- 



° Optimal estimation theory assumes that aU the errors are 
unbiased. 
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lution has not converged (determined by the method 
presented above within 10 iterations) then the calcu- 
lation is halted and a numerical 'flag' set. 

If, after convergence has been determined, the J(x) 
value is greater than the value given the number 
of degrees of freedom at a set confidence level (in this 
case 99.9%) then a numerical 'flag' is set. Retrievals 
with flags set are omitted from any statistics. 

Furthermore, to ensure the solution computed at 
each iteration is physical, the ln(specific humidity) el- 
ements are checked for super-saturation and corrected 
if necessary. 

4. Results 

In this section the performance of the optimal esti- 
mation retrieval scheme is examined using simulated 
profiles and realistic error estimates. For each simu- 
lated profile, a 'true' profile is established by taking 
one of a set of profiles of UKMO model analyses from 
which to compute the 'true' observation vector. The 
associated background atmospheric profile is calcu- 
lated by perturbing the 'true' profile thus 

n 

x'' = x*+^e,Ay'P, (12) 

1=1 

where the superscript t denotes the 'truth', Ai and 
are the ith eigenvalues and eigenvectors of the fore- 
cast error covariance matrix and represents the ith. 
number drawn from a normal distribution of random 
numbers. 

Observation noise is modelled and superimposed 
onto the observation vector using the method anal- 
ogous to equation ^ utilising the eigenvector and 
eigenvalues from the total observation error covari- 
ance matrix. 

4.1. Ensemble of numerical simulations 

Using the method described by equation |l^, sim- 
ulated observations and realistic background profiles 
were produced. These were used as inputs to the in- 
version scheme to obtain retrieved profiles which were 
subsequently compared with the 'true' profiles to as- 
sess the impact of the observations on the background 
information. 

The observation level values (i.e. impact parame- 
ter, local radius of curvature and geographical posi- 
tion) have been taken from data during ' prime-times 'Fl 



1 and 2 [Rocken et al, 1997], in an effort to simulate 
realistic retrieval scenarios. The latitudinal and longi- 
tudinal distribution of these occultation events (deter- 
mined by GPS sampling) are varied, thus providing a 
mixed ensemble of polar, tropical and mid-latitudinal 
occultation events. 

Five hundred profiles with random temperature, 
humidity and surface pressure conditions have been 
tested, and successful retrievals (i.e. which pass qual- 
ity control) were obtained in all but eight cases. In 
most cases convergence is obtained within three or 
four iterations. In general, the J(x) values at conver- 
gence were comparable to the number of degrees of 
freedom considered, as expected [Marks and Rodgers, 
1993]. This suggests that although the quantity is 
only strictly valid for linear problems it can be used 
reliably as a quality control for the retrieval. 

The eight cases which do not pass the specified 
quality control have been examined. They are found 
to be cases in which the cost function has been min- 
imised successfully but the converged value is too 
large compared with the distribution. These spu- 
rious converged profiles represent artificial outliers, 
which are generated when the increment described by 
equation |l^ is large enough to make the inverse prob- 
lem grossly non-linear. A small number of profiles are 
expected to have this problem due to the normal dis- 
tribution of random numbers used in the method of 
simulating atmospheric profiles. The profiles which 
failed the quality control have not been included in 
the statistics shown. 

For each successful retrieval, the retrieval error and 
the background error (first-guess error) have been cal- 
culated, and the mean and standard deviation values 
of these data have been computed. The standard de- 
viation values represent the errors ascribed to each 
element of the solution and background state vector, 
and can be compared directly with the square root 
values of the principal diagonal of the background er- 
ror covariance matrix assumed in the retrieval. The 
ratio of the retrieval error estimates to the forecast er- 
ror estimates is related to the amount of information 
the measurements supply to the NWP system. 

An improvement vector is defined which indicates 
how the retrieval has improved the knowledge of the 
background state throughout the atmospheric profile, 
and will be used to complement the r.m.s. statistics 
presented. The improvement vector is given by 

anti-spoofing military encryption. 
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J7,- = 100 X 



1 - 



(13) 



where j is the matrix and improvement vector element 
index, and all other variables are as before. 

The results from the ensemble of simulated re- 
trievals are summarised by figure HI The upper panels 
((a) and (b)) show the computed r.m.s. errors from 
the the simulated retrievals. The forecast errors re- 
semble those shown by figure || as expected, modified 
slightly by the modelled observation noise. 

The temperature improvement vector suggests that 
optimal estimation considerably improves upon the 
prior knowledge of the atmospheric temperature, from 
the lower-troposphere to the mid-stratosphere. The 
temperature improvement vector declines in the up- 
per stratosphere partly because of the upper noise 
limit used, which is comparable to forecast errors in 
observation space, and partly because it represents 
the lower limit of the observation vectors used in the 
ensemble of simulations. The temperature improve- 
ment vector declines in the lower-troposphere where 
humidity becomes more significant, and both quan- 
tities are retrieved simultaneously, thus the emphasis 
is shifted from temperature to ln(specific humidity). 
It is clear from the plot that there is a gradual de- 
cline in the temperature retrieval quality from 300 to 
1000 hPa as the humidity retrieval quality improves, 
where the background uncertainty information is be- 
ing used to resolve the temperature-specific humidity 
ambiguity in the refractivity. 

Both the temperature and humidity knowledge de- 
cline near the surface because the majority of occulta- 
tion events presented here terminate typically above 
one kilometre. 

The results shown by figure || confirm that the 
method is suitable for the purpose of non-linear op- 
timal estimation using RO measurements. Together, 
the theoretical r.m.s. errors and the computed im- 
provement vectors suggest that there are improve- 
ments in the prior knowledge of the atmospheric state 
from near-surface to the upper stratosphere. In par- 
ticular, the level of surface pressure improvement sug- 
gests that the RO observations can improve the prior 
knowledge of the surface pressure. 

4.2. Solution error characterisation 

Using the optimal estimation inverse theory out- 
lined in section |^ an error analysis can be obtained 



which allows a full characterisation of the solution (for 
a detailed account see Rodgers, [1990]). The error as- 
sociated with the solution vector can be split into its 
constituent parts, namely error from the background 
error estimates, forward modelling, forward model pa- 
rameters and observation noise. A mid-latitude re- 
trieval^ (which spans 0.7—60 km) has been used as 
an example to illustrate the method (figure ||). 

Panels (a) and (b) show that the a priori pro- 
vides almost all the information to the temperature 
and ln(specific humidity) retrieval above 10 hPa and 
500 hPa, respectively. 

Panels (c) and (d) show the observation noise con- 
tribution to the retrieval error which is comparatively 
small. 

Panels (e) and (f) show the forward modelling con- 
tribution to the retrieval error. The structure of this 
contribution is very similar to that of the observation 
noise, but has a larger associated error. In general, 
forward modelling represents the second largest con- 
tribution to retrieval error. 

The temperature solution error contributions shown 
by panels (c) and (e) peak at 3 hPa, at the point where 
the stratospheric noise limit contribution to the total 
observation error budget peaks (figure ||). The ob- 
servation contribution to the solution error begins to 
increase above about 0.3 hPa. Above this pressure 
level the a priori increases more rapidly, and so the 
observation is given more weight, resulting in a larger 
contribution to the solution error. 

Panels (g) and (h) show that forward model pa- 
rameter error represents the smallest contribution to 
the total retrieval error. This contribution to the tem- 
perature and ln(specific humidity) solutions peak near 
1000 hPa due to the local radius of curvature error. 
The contribution to the surface pressure solution rep- 
resents a small fraction of the total retrieval error. 

Figure ^ indicates that the dominant contributions 
to the solution error are from the a priori and forward 
modelling, suggesting that particular efforts should be 
made to improve their accuracies. It should be noted 
that the solution error depends on the assumptions 
made about the uncertainty statistics. 

4.3. Quality of surface pressure retrievals 

It is found that the quality of the surface pres- 
sure retrievals is dependent on the vertical extent of 



° Henceforth will be referred to as the example occultation 
profile. 
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the occultation, i.e. how closely it approaches the sur- 
face. To illustrate this point, the example occultation 
profile is used. By systematically removing observa- 
tions near the surface and re-retrieving temperature, 
ln(spccific humidity) and surface pressure, it can be 
shown how the vertical range of the occultation is im- 
portant to the quality of the surface pressure retrieval 
(figure |). 

Panels (a), (b) and (c) show the temperature, 
ln(specific humidity) and surface pressure improve- 
ment vectors get smaller in the lower atmosphere with 
increasing values for the lowest geometric height level, 
as expected. Panel (d) shows that the retrieval which 
includes all the observations provides enough informa- 
tion to retrieve accurately the true value for surface 
pressure; as the number of near-surface observations 
decreases the retrieval becomes smaller and smaller, 
approaching the prior value. 

It is interesting to note that the surface pressure 
information does not decrease as quickly as expected. 
Indeed there is still a considerable amount of surface 
pressure information towards the upper troposphere. 
This variation in surface pressure information is due 
to the link between height and pressure through the 
hydrostatic relation. 

It can be concluded from this experiment that us- 
ing optimal estimation, radio occultation measure- 
ments of the Earth possess surface pressure informa- 
tion even if the occultation has been completed in the 
mid-troposphere (e.g. due to atmospheric multipath 
interrupting the transmitted signal). This is impor- 
tant to appreciate when validating surface pressure 
retrievals using real data. 

4.4. Sensitivity to inverse model statistics 

In this section we present the results from a sensi- 
tivity study in which the statistics used to compute 
the optimal estimate are changed in order to investi- 
gate the retrieval sensitivity to such alterations. The 
three statistics that are altered are the observation 
noise, the forward modelling error and the forecast 
error since these represent the largest contributions 
to the solution error as shown by figure ^. 

Altering the observation noise may simulate the 
possible changes in observation noise sources, e.g. im- 
proved high-order ionospheric or poor quality clocks 
aboard the GPS receivers. 

Changing the forward modelling error is a crude 
method of simulating the possibility of assimilating 
intrinsically high resolution GPS radio occultation 



measurements with lower or higher resolution NWP 
model fields, and/or an occultation through a frontal 
system or a relatively horizontally inhomogeneous at- 
mosphere. 

Altering the prior error covariance matrix repre- 
sents the effect of changing the prior knowledge of 
the atmosphere. Increasing these error estimates can 
represent the extent of the knowledge of the atmo- 
sphere in the southern hemisphere where other at- 
mospheric information is sparse. Decreasing this er- 
ror may represent a more realistic model dynam- 
ics/climatology and/or greater confidence in other 
similar atmospheric observations used to initialise the 
model. For this particular test the variance values are 
changed, whilst retaining the existing correlations. 

Figure ^ shows the improvement vectors from the 
retrieval sensitivities outlined above. Improvement 
vectors are used to present the results of this study 
because it is the relative improvement on the prior es- 
timate of the atmospheric state that we are interested 
in. 

Panels (a) and (b) shows that doubling or halv- 
ing the observation noise has little effect on the over- 
all improvement, reflecting the contribution from this 
error source to the total observation error covariance 
matrix (figure ^). The changes due to ln(specific hu- 
midity) are very slight. The surface pressure improve- 
ment changes typically by a few percent. 

Panels (c) and (d) correspond to increasing or de- 
creasing the forward modelling error (by 50%). This 
error source provides the largest contribution to the 
total observation error, and as such has a large influ- 
ence on the degree of improvement. The difference 
in the temperature improvement is of the order of 
15% in the upper troposphere and lower stratosphere, 
above which other error sources are more important, 
and below which (the lower troposphere) the empha- 
sis is shifted from the temperature retrieval to the 
ln(specific humidity) and surface pressure retrieval. 
The improvement response for ln(specific humidity) 
is less pronounced than for temperature peaking at 
±8%. The surface pressure improvement variation is 
approximately ±15%. 

Panels (e) and (f) shows that increasing or de- 
creasing the standard deviation of the prior error in- 
creases or decreases the improvement of temperature, 
ln(specific humidity) throughout the vertical range of 
the observations as expected. Decreasing the a pri- 
ori error means better background knowledge, con- 
sequently the weighting of the a priori/observation 
is increased/decreased. This corresponds to a small 
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improvement relative to the background atmospheric 
knowledge. Positive and negative temperature im- 
provement differences are of the order of 20% through- 
out the range described by the observation vector; the 
ln(spccific humidity) improvement is of the order of 
±10%; and the surface pressure improvement is of the 
order of ±10%. 

This sensitivity study has looked at some of the ex- 
treme case scenarios in which the statistics assumed 
for the optimal estimation inverse method have been 
changed. It has been shown that the retrieval method 
is most sensitive to the background errors and forward 
modelling errors; the latter being related to errors due 
to the spherical symmetry assumption. The results 
from increasing the background error estimates are es- 
pecially important to note since they represent a real 
possibility when dealing with any reasonable measure- 
ments in the data sparse southern hemisphere. 

5. Conclusions 

In this paper we have demonstrated a prototype 
optimal estimation inverse method for GPS radio oc- 
cultation observations. 

The method is a non-linear, maximum a posteri- 
ori technique which can accommodate most aspects 
of the real radio occultation problem. In particular, it 
is able to account for some of the error incurred from 
assuming local spherical symmetry which is not pos- 
sible using the analytic inverse method. The optimal 
estimation technique handles the temperature— water 
vapour ambiguity in a more rigorous way, rather than 
the sub— optimal manner inherent with the analytic 
inverse method. 

The optimal estimation inverse method is used here 
as an iterative method but is found to be stable and 
to converge rapidly in most cases. The value of the 
cost function at each iteration can be used reliably 
to judge convergence and as an indicator of sensi- 
ble results, allowing anomalous retrieval conditions 
to be recognised and omitted, thus preventing gross 
retrieval errors. 

The method is shown to be suitable for retrieving 
values for surface pressure. Hence, this method of 
utilising radio occultation observations of the Earth's 
atmosphere has the potential to improve both atmo- 
spheric and oceanographic models, which may lead 
to improved predictions of the weather and climate. 
The quality of the surface pressure retrieval is shown 
to depend on the vertical extent of the occultation, 
i.e. higher quality retrievals are attainable with oc- 



cultations that reach low altitudes. 

It should be noted that the background statistics 
assumed in the paper represent global statistics. How- 
ever, for purposes of demonstrating a prototype re- 
trieval scheme for radio occultation observations they 
are found to be adequate. It has been shown that 
the retrieval accuracies, and hence the weight which 
should be given to the data in the subsequent model 
analysis, are sensitive to both forecast error uncer- 
tainties and values used to describe the forward mod- 
elling error. 
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Ray tangent 
point 




Figure 1. Defining the radio occultation geometry used to obtain bending angle information from the time delay 
caused by the Earth's atmosphere. A tangential sphere is superimposed on to the oblate Earth (exaggerated) to 
emphasise the position of the local radius of curvature at the ray pcriapsis. The dashed lines indicate motion. 
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Temperature uncertainty [K] ln(specifiG humidity) uncertainty 

Figure 2. The squaro-root values of the principal diagonal of the UKMO forecast error covariance matrix. Panels 
(a) and (b) show the assumed uncertainty for prior temperature and ln(spcc;ific humidity), respectively, with the 
assumed prior uncertainty for the surface pressure element inset of panel (b). 
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Figure 3. Measurement noise budget and total measurement error budget for the RO optimal estimation inverse 
method [after Luntama, [1997]]. The error estimates shown in panel (a) represent normal atmospheric conditions, 
i.e. small multipath error (3 mm) and normal ionospheric conditions. Panel (b) shows the percentage contributions 
from the different error sources to the total bending angle error. Panel (c) shows typical standard deviation values 
from the principal diagonal from each contribution to the total measurement error covariance matrix; and panel 
(d) shows the individual measurement error sources as a percentage proportion of the total error. 
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Figure 4. Theoretical error estimates from the ensemble of simulated retrievals. The upper panels show the 
theoretical r.m.s errors, where the solid and dotted lines represent the retrieval and background errors, respectively. 
Panel (a) shows the errors from the temperature solution and the panel (b) shows the errors from the In (specific 
humidity) solution with the surface pressure solution error inset. The lower panel show the corresponding improve- 
ment vector. Panels (c) and (d) show the temperature and ln(specific humidity) elements with the surface pressure 
improvement inset of panel (d). 
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Figure 5. Solution error characterisation of the optimal estimation inverse method for the example occultation 
profile. Left panels of a pair show the contribution to the temperature retrieval error and the right panels show 
the error contributions to the ln(specific humidity) retrieval error. The surface pressure retrieval error is inset of 
the left panels. Panels (a) and (b) show the a priori contribution to the retrieval error, where the dashed shown 
represent the standard deviation values of the principal diagonal of the a priori error covariance matrix; (c) and (d) 
the measurement noise contribution; (e) and (f) the forward modelling contribution; and (g) and (h) the forward 
model parameter contribution. The total solution error is superimposed on all panels (dotted line) to provide an 
indication as to which contributions are prominent. 
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Figure 6. Investigating the relationship between the quaHty of surface pressure retrievals and the lowest geometric 
height level of the occultation using the example occultation profile. Panel (a), (b) and (c) show the temperature, 
ln(specific humidity) and surface pressure improvement vectors; and panel (d) shows the difference (with errors 
bars) between the truth and the retrieval. The broken line in panel (d) represents the difference between the 'true' 
and the assumed value for the surface pressure. 




Figure 7. Improvement vectors from the retrieval sensitivity studies using the example occultation profile. In 
general, left panels show the temperature improvement vector elements and the right panels show the ln(specific 
humidity) improvement vector elements, with the surface pressure improvement vector clement inset. For panels 
(a) — (f), dashed curves represent a reduction in the quantity altered, dotted curves represent an increase in the quan- 
tity altered and solid curves represent the control case. Panels (a) and (b) show the results from doubling/halving 
the measurement noise; panels (c) and (d) show the results from increasing and decreasing the representative error 
by ±50%; and panels (e) and (f) show the results from increasing and decreasing the principal diagonal standard 
deviation errors from the a priori error covariance matrix by ±50%. 



